!/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!*  Copyleft (!c) 2004-2005 Div. Nucl. Med., Cyric, Tohoku Univ.    *
!*  Copyleft (!c) 2005-2011 Elseif Laboratory.                      *
!*  elseifkk@users.sf.net                                           *
!*                                                                  *
!*  All Rights Reversed.                                            *
!*                                                                  *
!*                                                                  *
!*  LICENSE AGREEMENT:                                              *
!*                                                                  *
!*  You agree that we don't agree to anything.                      *
!*                                                                  *
!*                                                                  *
!*  DISCLAIMER OF WARRANTIES:                                       *
!*                                                                  *
!*  This software is provided to you "AS IS," and the authors and   *
!*  contributors disclaim any warranty or liability obligations to  *
!*  you of any kind, implied or express, including but not limited  *
!*  to the warranties of merchantability, fitness for a particular  *
!*  purpose, and non-infringement.                                  *
!* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
!
! bitmap image utility; Daddy version.
!
! TTD:
!   none.
! note:
!   none.
!
module butyD
  implicit none
  real*8,parameter::PID=&
       & 6.283185307179586476925286766559005768394338798750211641949889184615
  real*8,parameter::PI=&
       & 3.141592653589793238462643383279502884197169399375105820974944592307
  real*8,parameter::PIH=&
       & 1.570796326794896619231321691639751442098584699687552910487472296153
  real*8,parameter::PIQ=&
       & .7853981633974483096156608458198757210492923498437764552437361480769
  real*8,parameter::SQRT2=&
       & 1.414213562373095048801688724209698078569671875376948073176679737990
  real*8,parameter::PVAL_OUT_OF_SIGHT = -TINY(0.0)
  real*8,parameter::PVAL_ALLOWABLE_NEGATIVE = -1000.
  integer*4,parameter::GMWRF_USE_SMOOTH_PDF = 1
  integer*4,parameter::GL_MAX = 256
  integer*4,parameter::GL_MIN = 63
  integer*4,parameter::MV_INT_TYPE_NN    =  0
  integer*4,parameter::MV_INT_TYPE_L2p   =  1
  !integer*4,parameter::MV_INT_TYPE_Lnp  =  n > 1
  integer*4,parameter::HIST_INT_TYPE_PV  = -3
  integer*4,parameter::MV_INT_TYPE_SINC  = -2
  integer*4,parameter::GMWRF_SMOOTH_PDF  = 1
  integer*4,parameter::DT_RKIND_4 = 32
  integer*4,parameter::DT_RKIND_8 = 64
  integer*4,parameter::DT_IKIND_1 = 2
  integer*4,parameter::DT_IKIND_2 = 4
  integer*4,parameter::DT_IKIND_4 = 8
  integer*4,parameter::DT_UNKNOWN = 0
  integer*4,parameter::DTF_RLC    = Z'01FF' ! run-length comp'ed
  integer*4,parameter::DT_ALL     = Z'00FF'
  integer*4,parameter::IMG_ALLOCED_R8 = Z'0001'
  integer*4,parameter::IMG_ALLOCED_R4 = Z'0002'
  integer*4,parameter::IMG_ALLOCED_I4 = Z'0004'
  integer*4,parameter::IMG_ALLOCED_I2 = Z'0008'
  integer*4,parameter::IMG_ALLOCED_I1 = Z'0010'
  integer*4,parameter::IMG_ALLOCED    = Z'001F'
  type image_data 
     integer*4 n(3)         ! number of pixels
     real*8 d(3)            ! pixel size
     real*8 xmin(3),xmax(3) ! image size 
     real*8 offset(3)       ! index to coord trans. 
     real*8 pmin,pmax       ! pixel value
     real*8 pave,pvar       ! pixel value average and var
     real*8 plo,pup         ! threshold for pixel value
     real*8 a,b             ! pixel value to index trans. i=a*p+b
     integer*4 ns           ! number of color index
     integer*4 stat         ! img stat
     integer*4 datatype     ! 4 or 8 (real4 or real8)
     real*8 cf              ! calibration factor
     real*8 H               ! entropy
     real*8,pointer::img_r8(:,:,:)
     real*4,pointer::img_r4(:,:,:)
     integer*1,pointer::img_i1(:,:,:)
     integer*2,pointer::img_i2(:,:,:)
     integer*4,pointer::img_i4(:,:,:)
  end type image_data
  !
  type Sinc_Table
     real*8,pointer::kt(:) ! kernel table
     integer*4 m(3)        ! mask size
     integer*4 n           ! allocate(kt(0:n))
     integer*4 d           ! sin(x)/x=kt(x*d)
  end type Sinc_Table
  !
  type LnP_Table
     integer*4 nh(3)
     integer*4 lo(3)
     integer*4 tlo(3)
     integer*4 tup(3)
     integer*4 kmax,kmin
     real*8,pointer::nu(:,:)
  end type Lnp_Table
  !
  type Interpolation_data
     !   integer*4 stat
     integer*4 it            ! interpolation type
     type(Sinc_Table) sinct  ! sinc int. table and all
     type(Lnp_Table) lnpt    ! Lagrange n point int. table and all
  end type Interpolation_data
  !
  type Gmwr_data
     integer*4 flag     ! controls gmwr function
     real*8 pdf_smtw    ! pdf gaussian smoothing width
     type(Interpolation_data) ipd
  end type Gmwr_data
  !
  type rigid_transform_param
     ! the tlansformed coordinate x' will be
     ! x'= MSR(x-d), where R, S and M represent rotation,
     ! sheer and magnification matrix and d is tlans. vec.
     real*8 d(3) ! tlanslation vector
     real*8 t(3) ! Euler rotation angle;
     real*8 m(3) ! magnification factor
     real*8 s(3) ! sheer factor
     ! actual rotation order is t2(z) -> t1(x) -> t3(z)
  end type rigid_transform_param
  !
  type smooth2d_table
     integer*4 n
     real*8,pointer::w(:,:)
  end type smooth2d_table
  !
  ! these are only meant as placeholder used by buty 
  ! routines; they must not be directly accessed.
  type(smooth2d_table),private,target::smt2dt
  type(Gmwr_data),pointer::bws_gmwrd
  type(Interpolation_Data),pointer::bws_ipd
  type(Sinc_Table),pointer::bws_sinct
  type(Lnp_Table),pointer::bws_lnpt
  real*8,pointer::bws_h(:,:,:) !work space for buty routines
  !------------------------------------------------------!
  interface operator(.by.)
     module procedure mv_coord
  end interface
  !------------------------------------------------------!
  interface operator(.invby.)
     module procedure mv_coord_inv
  end interface
  !------------------------------------------------------!
contains
  !------------------------------------------------------!
  subroutine uninit_gmwr(gmwrd)
    type(gmwr_Data),intent(inout),target::gmwrd
    call uninit_ipd(gmwrd%ipd)
    call uninit_smooth2d
    gmwrd%flag=0
  end subroutine uninit_gmwr
  !------------------------------------------------------!
  subroutine init_gmwr(gmwrd,int_type,&
       nx_lnpt,&
       ny_lnpt,&
       nz_lnpt,&
       nx_sinct,&
       ny_sinct,&
       nz_sinct,&
       d_sinct,&
       pdf_smtw,&
       pdf_smtw_th)
    type(gmwr_Data),intent(inout),target::gmwrd
    integer*4,intent(in)::int_type
    integer*4,intent(in),optional::nx_lnpt,ny_lnpt,nz_lnpt
    integer*4,intent(in),optional::nx_sinct,ny_sinct,nz_sinct,d_sinct
    real*8,intent(in),optional::pdf_smtw,pdf_smtw_th
    bws_gmwrd=>gmwrd
    gmwrd%flag=0
    call init_ipd(gmwrd%ipd,int_type,&
         nx_lnpt,&
         ny_lnpt,&
         nz_lnpt,&
         nx_sinct,&
         ny_sinct,&
         nz_sinct,&
         d_sinct)
    if(present(pdf_smtw).and.pdf_smtw/=0) then
       call init_smooth2d(pdf_smtw,pdf_smtw_th)
       gmwrd%pdf_smtw=pdf_smtw
       gmwrd%flag=ior(gmwrd%flag, GMWRF_USE_SMOOTH_PDF)
    endif
  end subroutine init_gmwr
  !
  subroutine uninit_ipd(ipd)
    type(Interpolation_Data),intent(inout),target::ipd
    select case(ipd%it)
    case(MV_INT_TYPE_SINC)
       call set_sinct(ipd%sinct)
       call regist_sinct
    case(MV_INT_TYPE_L2P,MV_INT_TYPE_NN,HIST_INT_TYPE_PV)
    case default
       call set_lnpt(ipd%lnpt)
       call regist_lnpt
    end select
  end subroutine uninit_ipd
  !
  subroutine init_ipd(ipd,int_type,&
       nx_lnpt,&
       ny_lnpt,&
       nz_lnpt,&
       nx_sinct,&
       ny_sinct,&
       nz_sinct,&
       d_sinct)
    type(Interpolation_Data),intent(inout),target::ipd
    integer*4,intent(in)::int_type
    integer*4,intent(in),optional::nx_lnpt,ny_lnpt,nz_lnpt
    integer*4,intent(in),optional::nx_sinct,ny_sinct,nz_sinct,d_sinct
    bws_ipd=>ipd
    ipd%it=int_type
    select case(int_type)
    case(MV_INT_TYPE_NN,MV_INT_TYPE_L2P)
       return
    case(MV_INT_TYPE_SINC)
       call set_sinct(ipd%sinct,nx_sinct,ny_sinct,nz_sinct,d_sinct)
       call regist_sinct(ipd%sinct)
    case(HIST_INT_TYPE_PV)
       ! ignore
    case default
       call set_lnpt(ipd%lnpt,nx_lnpt,ny_lnpt,nz_lnpt)
       call regist_lnpt(ipd%lnpt)
    end select
  end subroutine init_ipd
  !
  subroutine regist_lnpt(lnpt)
    type(LnP_Table),intent(in),target,optional::lnpt
    if(present(lnpt)) then
       bws_lnpt=>lnpt
    else
       nullify(bws_lnpt)
    end if
  end subroutine regist_lnpt
  !
  subroutine set_lnpt(lnpt,nx,ny,nz)
    type(LnP_Table),intent(inout)::lnpt
    integer*4,intent(in),optional::nx
    integer*4,intent(in),optional::ny,nz
    integer*4 i,k
    integer*4 n(3)
    if(.not.present(nx)) then
       call unset_lnpt
       return
    endif
    n=nx
    if(present(ny)) n(2)=ny
    if(present(nz)) n(3)=nz
    where(iand(n,1)==0)
       lnpt%nh=n/2
       lnpt%lo=-lnpt%nh+1
       lnpt%tlo=1
       lnpt%tup=n
    elsewhere
       lnpt%nh=(n-1)/2
       lnpt%lo=-lnpt%nh
       lnpt%tlo=0
       lnpt%tup=n-1
    endwhere
    lnpt%kmax=max(lnpt%nh(1),lnpt%nh(2),lnpt%nh(3))
    lnpt%kmin=min(lnpt%lo(1),lnpt%lo(2),lnpt%lo(3))
    allocate(lnpt%nu(lnpt%kmin:lnpt%kmax,3))
    do i=1,3
       do k=lnpt%kmin,lnpt%kmax
          lnpt%nu(k,i)=1/(fact(-lnpt%lo(i)+k)*fact(lnpt%nh(i)-k))
          if(iand(lnpt%nh(i)+k,1)/=0) lnpt%nu(k,i)=-lnpt%nu(k,i)
       end do
    end do
  contains
    real*8 function fact(xx)
      integer*4,intent(in)::xx
      integer*4 ii
      fact=1
      if(xx/=0) then
         do ii=1,xx
            fact=fact*ii
         end do
      end if
    end function fact
    !
    subroutine unset_lnpt
      lnpt%kmax=0
      lnpt%kmin=0
      if(associated(lnpt%nu)) deallocate(lnpt%nu)
    end subroutine unset_lnpt
  end subroutine set_lnpt
  !
  subroutine regist_buty_work_space(h)
    real*8,intent(in),target,optional::h(:,:,:)
    if(present(h)) then
       bws_h=>h
    else
       nullify(bws_h)
    end if
  end subroutine regist_buty_work_space
  !
  subroutine regist_sinct(sinct)
    type(Sinc_Table),intent(in),target,optional::sinct
    if(present(sinct)) then
       bws_sinct=>sinct
    else
       nullify(bws_sinct)
    end if
  end subroutine regist_sinct
  !
  subroutine set_sinct(sinct,nx,ny,nz,d)
    type(Sinc_Table),intent(inout)::sinct
    integer*4,intent(in),optional::nx
    integer*4,intent(in),optional::ny,nz
    integer*4,intent(in),optional::d
    real*8 dx
    integer*4 i,nmax
    if(.not.present(nx)) then
       call unset_sinct
       return
    end if
    sinct%m=nx
    if(present(ny)) sinct%m(2)=ny
    if(present(nz)) sinct%m(3)=nz
    nmax=max(sinct%m(1),sinct%m(2),sinct%m(3))
    sinct%n=(1+nmax)*d
    sinct%d=d
    allocate(sinct%kt(0:sinct%n))
    if(present(d)) then
       dx=1/dble(d)
    else
       dx=1
    endif
    do i=0,sinct%n
       sinct%kt(i)=sinc_kernel(dble(i)*dx)
    end do
  contains
    real*8 function sinc_kernel(xin)
      real*8,intent(in)::xin
      if(xin/=0) then
         sinc_kernel=dsin(PI*xin)/(PI*xin)
      else
         sinc_kernel=1
      endif
    end function sinc_kernel
    !
    subroutine unset_sinct
      if(associated(sinct%kt)) deallocate(sinct%kt)
      sinct%n=0
      sinct%d=0
    end subroutine unset_sinct
  end subroutine set_sinct
  !======================================================!  
  function mv_coord(a,x)
    ! a must be reduced affine matrix
    real*8 mv_coord(3)
    real*8,intent(in)::a(3,4)
    real*8,intent(in)::x(3)
    mv_coord=(x(1)-a(1,4))*a(:,1)&
         +(x(2)-a(2,4))*a(:,2)&
         +(x(3)-a(3,4))*a(:,3)
  end function mv_coord
  function mv_coord_inv(a,x)
    ! a must be reduced affine matrix
    real*8 mv_coord_inv(3)
    real*8,intent(in)::a(3,4)
    real*8,intent(in)::x(3)
    mv_coord_inv=x(1)*a(1,1:3)&
         +x(2)*a(2,1:3)&
         +x(3)*a(3,1:3)+a(:,4)
  end function mv_coord_inv
  !------------------------------------------------------!
  subroutine affine_rigid_to_rtp(a,rtp)
    real*8,intent(out)::a(3,4)
    type(rigid_transform_param),intent(out)::rtp
    real*8 c1,s1
    rtp%d(1:3)=a(1:3,4)
    c1=a(3,3)
    if(abs(c1)/=1) then
       rtp%t(3)=atan(a(1,3)/a(2,3))
       s1=dsqrt(1-a(3,3)**2)
       rtp%t(1)=asin(s1)
       rtp%t(2)=acos(a(3,2)/s1)
    else ! s1=0, |c1|=1, a23=a13=a31=0
       if(a(1,2)==-a(2,2).and.a(1,1)==-a(2,1))then
          !      c3=0 ! |s3|=1
          rtp%t(2)=atan((a(1,2)-a(2,2))/(a(1,1)-a(2,1)))
          rtp%t(3)=pih
       elseif(a(1,2)==-a(2,2)) then
          !      c2=0 ! |s2|=1
          rtp%t(2)=pih
          rtp%t(3)=-asin(a(2,2))
       elseif(a(1,1)==-a(2,1)) then
          !      s2=0 ! |c2|=1
          rtp%t(2)=0   
          rtp%t(3)=acos(a(1,1))
       else
          rtp%t(2)=-atan((a(1,1)+a(2,1))/(a(1,2)+a(2,2)))
          rtp%t(3)=acos(dsqrt((a(1,2)-a(2,2))**2+(a(1,1)-a(2,1))**2)/sqrt2)+piq
       endif
    endif
  end subroutine affine_rigid_to_rtp
  !------------------------------------------------------!
  subroutine set_affine_rigid(rtp,a)
    type(rigid_transform_param),intent(in)::rtp
    real*8,intent(out)::a(3,4)
    real*8 c(3),s(3)
    a(1:3,4)=rtp%d(1:3)
    c=cos(rtp%t)
    s=sin(rtp%t)
    a(1,1)=(c(2)*c(3)-c(1)*s(2)*s(3))
    a(2,1)=(-c(2)*s(3)-c(1)*s(2)*c(3))
    a(3,1)=(s(1)*s(2))
    a(1,2)=(s(2)*c(3)+c(1)*c(2)*s(3))
    a(2,2)=(-s(2)*s(3)+c(1)*c(2)*c(3))
    a(3,2)=(-s(1)*c(2))
    a(1,3)=(s(1)*s(3))
    a(2,3)=(s(1)*c(3))
    a(3,3)=(c(1))
  end subroutine set_affine_rigid
  !------------------------------------------------------!
  subroutine set_affine(rtp,a)
    type(rigid_transform_param),intent(in)::rtp
    real*8,intent(out)::a(3,4)
    call set_affine_rigid(rtp,a)
    ! magnification first
    a(1,:)=rtp%m(1)*a(1,:)
    a(2,:)=rtp%m(2)*a(2,:)
    a(3,:)=rtp%m(3)*a(3,:)
    ! then apply sheer
    a(1,:)=a(1,:)+rtp%s(1)*a(2,:)+rtp%s(2)*a(3,:)
    a(2,:)=a(2,:)+rtp%s(3)*a(3,:)
  end subroutine set_affine
  !------------------------------------------------------!
  subroutine print_image_data(x,title)
    use cuty
    character*128 str
    type(image_data),intent(in)::x
    character*(*),intent(in),optional::title
    integer*4 l
    if(present(title)) then
       l=conc(str,conc(str,conc(str,1,"=== "),title)," ===")
       call messn(str(1:l-1))
    endif
    !
    l=conc(str,conc(str,conc(str,conc(str,conc(str,1,&
         & "Matrix kind: "),x%datatype)," ("),trim(dt_to_str(x%datatype))),")")
    call messn(str(1:l-1))
    !
    l=conc(str,conc(str,conc(str,conc(str,conc(str,conc(str,1,&
         & "Dim: "),x%n(1)),"x"),x%n(2)),"x"),x%n(3))
    call messn(str(1:l-1))
    !
    l=conc(str,conc(str,conc(str,conc(str,conc(str,conc(str,1,&
         & "Vox. size: "),x%d(1)),"x"),x%d(2)),"x"),x%d(3))
    call messn(str(1:l-1))
    !
    l=conc(str,conc(str,conc(str,conc(str,conc(str,conc(str,conc(str,conc(str,&
         & conc(str,conc(str,conc(str,conc(str,1,&
         & "Img. size: "),x%xmin(1)),"->"),x%xmax(1)),", "),x%xmin(2)),"->"),&
         & x%xmax(2)),", "),x%xmin(3)),"->"),x%xmax(3))
    call messn(str(1:l-1))
    !
    l=conc(str,conc(str,1,"CF: "),x%cf)
    call messn(str(1:l-1))
    !
    l=conc(str,conc(str,conc(str,conc(str,conc(str,conc(str,conc(str,conc(str,1,&
         & "Cnt. min/max/ave/var: "),&
         & x%pmin),"/"),x%pmax),"/"),x%pave),"/"),x%pvar)
    call messn(str(1:l-1))  
    !
    l=conc(str,conc(str,conc(str,conc(str,1,&
         & "Th. lower/upper: "),x%plo),"/"),x%pup)
    call messn(str(1:l-1))
    !
    l=conc(str,conc(str,1,"ns: "),x%ns)
    call messn(str(1:l-1))
    !
  end subroutine print_image_data
  !------------------------------------------------------!
  subroutine init_image_data(x)
    type(image_data),intent(inout)::x
    call dealloc_image_data(x)
    x%stat=0
    x%n=0
    x%d=0
    x%xmin=0
    x%xmax=0
    x%offset=0
    x%pmin=0
    x%pmax=0
    x%pave=0
    x%pvar=0
    x%plo=-HUGE(0.d0)
    x%pup=HUGE(0.d0)
    x%a=0
    x%b=0
    x%ns=0
    X%cf=0
    x%datatype=0
  end subroutine init_image_data
  !------------------------------------------------------!
  subroutine cp_image_data(a,b)
    type(image_data),intent(in)::a
    type(image_data),intent(out)::b
    b%n=a%n
    b%d=a%d
    b%xmin=a%xmin
    b%xmax=a%xmax
    b%offset=a%offset
    b%pmin=a%pmin
    b%pmax=a%pmax
    b%pave=a%pave
    b%plo=a%plo
    b%pup=a%pup
    b%a=a%a
    b%b=a%b
    b%ns=a%ns
    b%cf=a%cf
    b%datatype=a%datatype
  end subroutine cp_image_data
  !------------------------------------------------------!
  subroutine dealloc_image_data(x)
    type(image_data),intent(inout)::x  
    if(iand(x%stat,IMG_ALLOCED_R8)/=0) deallocate(x%img_r8)
    if(iand(x%stat,IMG_ALLOCED_R4)/=0) deallocate(x%img_r4) 
    if(iand(x%stat,IMG_ALLOCED_I4)/=0) deallocate(x%img_i4)
    if(iand(x%stat,IMG_ALLOCED_I2)/=0) deallocate(x%img_i2) 
    if(iand(x%stat,IMG_ALLOCED_I1)/=0) deallocate(x%img_i1)
    x%stat=iand(x%stat,not(ior(IMG_ALLOCED_R8,&
         ior(IMG_ALLOCED_R4,&
         ior(IMG_ALLOCED_I4,&
         ior(IMG_ALLOCED_I2,&
         IMG_ALLOCED_I1))))))
    x%datatype=0
  end subroutine dealloc_image_data
  !------------------------------------------------------!
  subroutine alloc_image_data(x,n,kind)
    type(image_data),intent(inout)::x
    integer*4,intent(in),optional::n(3)
    integer*4,intent(in),optional::kind
    integer*4 istat,k
    call dealloc_image_data(x)
    if(.not.present(n)) return
    if(present(kind))then
       k=kind
    else
       k=DT_RKIND_8
    endif
    select case(k)
    case(DT_RKIND_4)
       allocate(x%img_r4(n(1),n(2),n(3)),stat=istat)
       if(istat==0) x%stat=ior(x%stat,IMG_ALLOCED_R4)
    case(DT_RKIND_8)
       allocate(x%img_r8(n(1),n(2),n(3)),stat=istat)
       if(istat==0) x%stat=ior(x%stat,IMG_ALLOCED_R8)
    case(DT_IKIND_4)
       allocate(x%img_i4(n(1),n(2),n(3)),stat=istat)
       if(istat==0) x%stat=ior(x%stat,IMG_ALLOCED_I4)
    case(DT_IKIND_2)
       allocate(x%img_i2(n(1),n(2),n(3)),stat=istat)
       if(istat==0) x%stat=ior(x%stat,IMG_ALLOCED_I2)
    case(DT_IKIND_1) 
       allocate(x%img_i1(n(1),n(2),n(3)),stat=istat)
       if(istat==0) x%stat=ior(x%stat,IMG_ALLOCED_I1)
    end select
    if(istat==0) then
       x%datatype=k
       x%n=n
    else
       x%n=0
    endif
  end subroutine alloc_image_data
  !------------------------------------------------------!
  subroutine print_affine(a)
    use cuty
    real*8,intent(in)::a(3,4)
    character*128 str
    integer*4 l,i,j
    do j=1,3
       l=1
       do i=1,4
          l=l+f2str(a(j,i),str(l:),FP2A_DEFAULT)
          str(l:l)=" "
          l=l+1
       enddo
       str(l:l)=char(10)
       call mess(str(1:l))
    end do
  end subroutine print_affine
  !------------------------------------------------------!
  subroutine inverse_affine(a)
    real*8,intent(inout)::a(3,4)
    real*8 s(3)
    s(1:2)=a(2:3,1)
    s(3)=a(3,2)
    a(2:3,1)=a(1,2:3)
    a(3,2)=a(2,3)
    a(1,2:3)=s(1:2)
    a(2,3)=s(3)
    s=a(:,4)
    a(:,4)=-(a(:,1)*s(1)+a(:,2)*s(2)+a(:,3)*s(3))
  end subroutine inverse_affine
  !------------------------------------------------------!
  subroutine expand_affine(a)
    real*8,intent(inout)::a(3,4)
    real*8 s(3)
    s=a(:,4)
    a(:,4)=(a(:,1)*s(1)+a(:,2)*s(2)+a(:,3)*s(3))
  end subroutine expand_affine
  !------------------------------------------------------!
  subroutine get_calibration_factor(id,dt,cf,r)
    type(image_data),intent(in)::id
    integer*4,intent(in)::dt
    real*8,intent(out)::cf
    real*8,intent(in),optional::r
    real*8 pmin,pmax,cfl,cfr,rr
    if(present(r))then
       rr=r
    else
       rr=1
    endif
    select case(dt)
    case(DT_IKIND_1)
       pmin=0
       pmax=255
    case(DT_IKIND_2)
       pmin=-32767
       pmax= 32767
    case(DT_IKIND_4)
       pmin=-2147483647
       pmax= 2147483647
    case(DT_RKIND_4,DT_RKIND_8)
       cf=1
       return
    end select
    cfl=id%pmin/(pmin*r)
    cfr=id%pmax/(pmax*r)
    cf=max(abs(cfl),abs(cfr))
  end subroutine get_calibration_factor
  !------------------------------------------------------!
  subroutine uninit_smooth2d()
    if(associated(smt2dt%w)) deallocate(smt2dt%w)
    smt2dt%n=0
  end subroutine uninit_smooth2d
  !------------------------------------------------------!
  subroutine init_smooth2d(sigma,threshold)
    real*8,intent(in)::sigma
    real*8,intent(in),optional::threshold
    real*8 th
    integer*4 ix,iy
    call uninit_smooth2d
    if(present(threshold)) then 
       th=threshold
    else
       th=3.d0
    endif
    smt2dt%n=max(1,nint(th*sigma))
    allocate(smt2dt%w(-smt2dt%n:smt2dt%n,-smt2dt%n:smt2dt%n))
    do iy=-smt2dt%n,smt2dt%n
       do ix=-smt2dt%n,smt2dt%n
          smt2dt%w(ix,iy)=g(sqrt(dble(ix**2+iy**2)),sigma)
       end do
    end do
    smt2dt%w=smt2dt%w/sum(smt2dt%w)
  contains
    real*8 function g(x,sgm)
      real*8,intent(in)::x
      real*8,intent(in)::sgm
      g=1/(sqrt(2*pi)*sgm)*exp(-((x/sgm)**2)/2)
    end function g
  end subroutine init_smooth2d
  !------------------------------------------------------!
  subroutine smooth2d(nx,ny,in,out)
    ! for pdf smoothing where the integral must be preserved
    integer*4,intent(in)::nx,ny
    real*8,intent(in)::in(nx,ny)
    real*8,intent(out)::out(nx,ny)
    real*8,parameter::G_THRESHOLD = 3
    integer*4 ix,iy,jx,jy,iix,iiy
    real*8,pointer::w(:,:)
    integer*4,pointer::ns
    real*8 s,sin,sout,norm
    ns=>smt2dt%n
    w=>smt2dt%w
    sin=sum(in)
    sout=0
    do iy=1,ny
       do ix=1,nx
          s=0
          do iiy=-ns,ns
             jy=iy+iiy
             if(jy<=0) cycle
             if(jy>ny) exit
             do iix=-ns,ns
                jx=ix+iix
                if(jx<=0) cycle
                if(jx>nx) exit
                s=s+w(iix,iiy)*in(jx,jy)
             end do
          end do
          out(ix,iy)=s
          sout=sout+s
       end do
    end do
    if(sout/=0) then
       norm=sin/sout
       out=out*norm
    endif
  end subroutine smooth2d
  !------------------------------------------------------!
  integer*4 function find_bg_level(n,h,lo)
    integer*4,intent(in)::n
    real*8,intent(in)::h(n)
    real*8,intent(in)::lo ! threshold relative to hmax
    real*8 hmax,th
    integer*4 imax,i
    imax=0
    hmax=-HUGE(0.d0)
    do i=1,n
       if(h(i)>hmax) then
          imax=i
          hmax=h(i)
       end if
    end do
    th=(1-lo)*hmax
    find_bg_level=0
    do i=imax,n
       if(h(i)<th) then
          find_bg_level=i
          exit
       end if
    end do
    FIND_BG_LEVEL=10
  end function find_bg_level
  !------------------------------------------------------!
  real*8 function get_max_d(d1,d2)
    real*8,intent(in)::d1(3)
    real*8,intent(in),optional::d2(3)
    get_max_d=max(d1(1),d1(2),d1(3))
    if(present(d2)) get_max_d=max(get_max_d,d2(1),d2(2),d2(3))
  end function get_max_d
  !------------------------------------------------------!
  subroutine get_gmask_size(sigma,th,n)
    real*8,intent(in)::sigma(3)
    real*8,intent(in)::th
    integer*4,intent(out)::n(3)
    where(sigma/=0)
       n=max(1,nint(sigma*th/2+0.5))
    elsewhere
       n=0
    end where
  end subroutine get_gmask_size
  !------------------------------------------------------!
  subroutine get_gmask(sigma,n,mask)
    real*8,intent(in)::sigma(3)
    integer*4,intent(in)::n(3)
    real*8,intent(out)::mask(-n(1):n(1),-n(2):n(2),-n(3):n(3))
    integer*4 ix,iy,iz
    do iz=0,n(3)
       do iy=0,n(2)
          do ix=0,n(1)
             mask(ix,iy,iz)=g(dble(ix),dble(iy),dble(iz))
          end do
! miser code BEGIN
          do ix=-n(1),-1 
             mask(ix,iy,iz)=mask(-ix,iy,iz)
          end do
       end do
       do iy=-n(2),-1
          do ix=-n(1),n(1)
             mask(ix,iy,iz)=mask(ix,-iy,iz)
          end do
       end do
    end do
    do iz=-n(3),-1
       do iy=-n(2),n(2)
          do ix=-n(1),n(1)
             mask(ix,iy,iz)=mask(ix,iy,-iz)
          end do
       end do
    end do
! miser code END
    mask=mask/sum(mask)
  contains
    real*8 function g(x,y,z)
      real*8,intent(in)::x,y,z
      g=1
      if(sigma(1)/=0) g=g*exp(-(x/sigma(1))**2/2)
      if(sigma(2)/=0) g=g*exp(-(y/sigma(2))**2/2)
      if(sigma(3)/=0) g=g*exp(-(z/sigma(3))**2/2)
    end function g
  end subroutine get_gmask
  !
  character*8 function dt_to_str(i)
    integer*4,intent(in)::i
    select case(i)
    case(DT_IKIND_1)
       dt_to_str="int1"
    case(DT_IKIND_2)
       dt_to_str="int2"
    case(DT_IKIND_4)
       dt_to_str="int4"
    case(DT_RKIND_4)
       dt_to_str="real4"
    case(DT_RKIND_8)
       dt_to_str="real8"
    case default
       dt_to_str="???"
    end select
  end function dt_to_str
end module butyD
